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Abstract 



Least-squares optimized polynomials are discussed which are needed in the two- 
step multi-bosonic algorithm for Monte Carlo simulations of quantum field theories 
• with fermions. A recurrence scheme for the calculation of necessary coefficients in 



the recursion and for the evaluation of these polynomials is introduced. 

1 Introduction 

In popular Monte Carlo simulation algorithms for QCD and other similar quantum field 
theories the main difficulty is the evaluation of the determinant of the fermion action 
matrix. This can be achieved by stochastic procedures with the help of auxiliary bosonic 
"pseudofermion" fields. 

In the two-step multi-bosonic (TSMB) algorithm 0] an approximation of the fermion 
determinant is achieved by the pseudofermion fields corresponding to a polynomial ap- 
proximation of some negative power a;~° of the fermion matrix . The auxiliary bosonic 
fields are updated according to the multi-bosonic action P|. The error of the polynomial 

*Talk given at the workshop on Numerical Challenges in Lattice QCD, August 1999, Wuppertal 
University; to appear in the proceedings. 
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approximation is corrected in a global accept-reject decision by using better polynomial 
approximations. Sometimes a reweighting of gauge configurations in the evaluation of ex- 
pectation values is also necessary. This can also be performed by high order polynomials. 

The polynomials used in the TSMB algorithm have to approximate the function 
x~'^P{x) in some non- negative interval x G [e, A], < e < A. Here P{x) is a known 
polynomial, typically another cruder approximation of x^". The approximation scheme 
and optimization procedure can be chosen differently. The least-squares optimization 
is an efficient and flexible possibility. (For other approximation schemes see 

In this review the basic relations for least-squares optimized polynomials are presented 
as introduced in |^. Particular attention is paid to a recurrence scheme which can be 
applied for determining the necessary high order polynomials and for evaluating them 
numerically. The details of the TSMB algorithm will not be considered. For a compre- 
hensive summary and references see 0. For experience on the application of TSMB in a 
recent large scale numerical simulation see |^ . 

2 Basic relations 

Least-squares optimization provides a general and flexible framework for obtaining the 
necessary optimized polynomials in multi-bosonic fermion algorithms. Here we introduce 
the basic formulae in the way it has been done in 0, |^ . 

We want to approximate the real function f{x) in the interval x G [e, A] by a polynomial 
Pn{x) of degree n. The aim is to minimize the deviation norm 



Here w{x) is an arbitrary real weight function and the overall normalization factor N^^x 
can be chosen by convenience, for instance, as 



A typical example of functions to be approximated is f\x) = x~"/-P(x) with a > and 
some polynomial P{x). The interval is usually such that < e < A. For optimizing the 
relative deviation one takes a weight function w{x) = f{x)~^. 

5n is a quadratic form in the coefficients of the polynomial which can be straightfor- 
wardly minimized. Let us now consider, for simplicity, only the relative deviation from 
the simple function f{x) = x~°' = w(a;)~^. Let us denote the polynomial corresponding 
to the minimum of 5„ by 




(1) 




(2) 



n 




(3) 
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Performing the integral in 5^ term by term we obtain 



5l = l-2Y.c.V!;-^+ Y: c.M:l,c.., (4) 

v=0 v\,V2=0 



where 

\ 1 _1_/-V-1-'H 1J 



y{") 



(A - e)(l + a + n-p) 



\l+2a+2ri~u\~U2 A+2a+2ri~u\—U2 

{\-e){l + 2a + 2n-Vi-V2) ' ^ ' 

The coefficients of the polynomial corresponding to the minimum of 5^, or of 5^, are 

n 

c. ^ c„.(a; e, A) = 5: M(^- V,(") . (6) 

1/1=0 

The value at the minimum is 

n 

5l^5l{a-e^) = l- K^r^M^^Vj,") . (7) 

The solution of the quadratic optimization in (^-(^ gives in principle a simple way to 
find the required least-squares optimized polynomials. The practical problem is, however, 
that the matrix M^"^ is not well conditioned because it has eigenvalues with very different 
magnitudes. In order to illustrate this let us consider the special case (a = l, A = l, e = 
0) with n = 10. In this case the eigenvalues are: 

0.4435021205e - 14 , 0.1045947635e - 11 , 0.1143819915e - 9 , 

0.7698917100e - 8 , 0.3571195735e - 6 , 0.1211873623e - 4 , 

0.3120413130e - 3 , 0.6249495675e - 2 , 0.9849331094e - 1 , 
1.075807246 . (8) 

A numerical investigation shows that, in general, the ratio of maximal to minimal eigen- 
values is of the order of (9(10^'^"). It is obvious from the structure of M*^"-* in (^) that 
a rescaling of the interval [e. A] does not help. The large differences in magnitude of 
the eigenvalues implies through large differences of magnitude in the coefficients c„y 
and therefore the numerical evaluation of the optimal polynomial P„(x) for large n is 
non-trivial. 

Let us now return to the general case with arbitrary function f{x) and weight w{x). 
It is very useful to introduce orthogonal polynomials (/i = 0, 1, 2, . . .) satisfying 

dxw{x)'^^^{x)^u{x) = 5^yqu . (9) 



and expand the polynomial in terms of them: 

n 

Pn{x) = Y^dnuM^) ■ (10) 
1^=0 

Besides the normalization factor gj, let us also introduce, for later purposes, the integrals 
Pj, and by 

f-A 



rX 

Pj, = I dxw{xy^^{xyx , 



A 

dxw{xyx'^ . (11) 

It can be easily shown that the expansion coefficients dm, minimizing 5„ are indepen- 
dent of n and are given by 

dnu = di,^ — , (12) 

where 

rX 

K = J dxw{xff{x)<i>^{x). (13) 

The minimal value of is 

n 

Rescaling the variable x by x' = px allows for considering only standard intervals, say 
[e/A, 1]. The scaling properties of the optimized polynomials can be easily obtained from 
the definitions. Let us now again consider the simple function f{x) = x~°' and relative 
deviation with w{x) ~ x°' when the rescaling relations are: 

Pn{a;ep,\p;x) = p~'^Pn{a;e,X;x/p) , 

Cnu{a;ep,Xp) = p''~"~"c„,,(Q;;e, A) . (15) 

In applications to multi-bosonic algorithms for fermions the decomposition of the 
optimized polynomials as a product of root-factors is needed. This can be written as 

n 

Pn{a;€,X;x) ^ Cno{a;e,X)'[l[x -rnj{a;e,X)] . (16) 

The rescaling properties here are: 

Cno{a; ep, Xp) = c„o(q;; e, A) , 

rnj{a;ep,Xp) = prnj{a;e,X) . (17) 
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The root-factorized form (0) can also be used for the numerical evaluation of the poly- 
nomials with matrix arguments if a suitable optimization of the ordering of roots is per- 
formed 0. 

The above orthogonal polynomials satisfy three-term recurrence relations which are 
very useful for numerical evaluation. In fact, at large n the recursive evaluation of the 
polynomials is numerically more stable than the evaluation with root factors. For general 
f{x) and w{x), the first two ortogonal polynomials with /i = 0, 1 are given by 

$o(a;) = 1 , $i(x) = x-— . (18) 

■50 

The higher order polynomials $^(a;) for = 2, 3, . . . can be obtained from the recurrence 
relation 

= (x + /?^)$^(x)+7^_i$^_i(x) , (/i = 1,2,...) , (19) 
where the recurrence coefficients are given by 

= , 7.-1 = . (20) 

Defining the polynomial coefficients f^y (0 < < /i) by 

^,{^) = E U^^'-" (21) 

the above recurrence relations imply the normalization convention 

^0 = 1, (/i = 0,l,2,...) . (22) 

The rescaling relations for the orthogonal polynomials easily follow from the defini- 
tions. For the simple function f{x) = x~" and relative deviation with w{x) = we 
have 

$^(a; pe, pX; x) = ^^(a; e. A; x/p) . (23) 
For the quantities introduced in (|1T]) this implies 

g,(a;pe,pA) = p'"+'+''' q^ia; e, X) , 
pAc^;pe,pX) = p2"+2+2>.(a;e,A) , 

s,(a;pe,pA) = p2°+i+'^ s,(a; e. A) . (24) 
For the expansion coefficients defined in ([T2|)-(|T3|) one obtains 

6,(«;pe,pA) = p^+i^'^ 6,(«; e. A) , 

dy{a; pe, pX) = p'"''" dy{a; e, X) , (25) 



and the recurrence coefficients in (p^ -(pO|) satisfy 



7^_i(a;pe,pA) = 7^_i(a; e, A) . (26) 

For general intervals [e, A] and/or functions f{x) = x~'^P{x) the orthogonal polynomi- 
als and expansion coefficients have to be determined numerically. In some special cases, 
however, the polynomials can be related to some well know ones. An example is the 
weight factor 

w^P^''\xf = {x-tY{\-xY . (27) 

Taking, for instance, p = 2a, a = this weight is similar to the one for relative deviation 
from the function /(x) = which would be just x^". In fact, for e = these are 
exactly the same and for small e the difference is negligible. The corresponding orthogonal 
polynomials are simply related to the Jacobi polynomials namely 

(X) = (A - e) V! I(P±^±^p(..P) I ^^-^-^ \ (28) 

Comparing different approximations with different (p, cr) the best choice is usually p = 
2a, (7 = which corresponds to optimizing the relative deviation (see the appendix of 

i)- 

For large condition numbers A/e least-squares optimization is much better than the 
Chebyshev approximation used for the approximation of x~^ in The Chebyshev 
approximation is minimizing the maximum of the relative deviation 

R{x) = xP{x) - 1 . (29) 

For the deviation norm 

5max = max \R{x)\ (30) 

the least-squares approximation is slightly worse than the Chebyshev approximation. An 
example is shown by fig. [1|. In the left lower corner the Chebyshev approximation has 
i?c(0.0002) = -0.968 compared to i?o(0.0002) = -0.991 for the least-squares optimiza- 
tion. For smaller condition numbers the Chebyshev approximation is not as bad as is 
shown by fig. |l]. Nevertheless, in QCD simulations in sufficiently large volumes the con- 
dition number is of the order of the light quark mass squared in lattice units which can 
be as large as 0(10^ - 10^). 

Figure |I| also shows that the least-squares optimization is quite good in the minimax 
norm in (0), too. It can be proven that 

\Ro{e)\ = max \Ro{x)\ (31) 
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hence the minimax norm can also be easily obtained from 



5^1= max \Ro{x)\ = \Ro{e)\ . (32) 

Therefore the least squares-optimization is also well suited for controlling the minimax 
norm, if for some reason it is required. 

In QCD simulations the inverse power to be approximated [a) is related to the number 
of Dirac fermion flavours: a = Nf/2. If only u- and d-quarks are considered we have 
Nf = 2 and the function to be approximated is x~^. The dependence of the (squared) 
least-squares norm in (|I]) on the polynomial order n is shown by fig. |^ for different values 
of the condition number A/e. The dependence on a = Nf/2 is illustrated by fig. 

Another possible application of least-squares optimized polynomials is the numerical 
evaluation of the zero mass lattice action proposed by Neuberger H^]. If one takes, 



for instance, the weight factor in ( ^71) corresponding to the relative deviation, then the 
function x~^^'^ has to be expanded in the Jacobi polynomials p(^'°). 



3 Recurrence scheme 

The expansion in orthogonal polynomials is very useful because it allows for a numerically 
stable evaluation of the least-squares optimized polynomials by the recurrence relation 
([T9|) . The orthogonal polynomials themselves can also be determined recursively. 

A recurrence scheme for obtaining the recurrence coefficients /3^,7^_i and expansion 
coefficients d^, = K/qu has been given in 0- In order to obtain qu,Pu contained in ( pO]) 
one can use the relations 

fi fj, 

(33) 

u=0 u=Q 

The coefficients themselves can be calculated from /n = —si/sq and (|1^) which gives 



f (1+1,1 = fti,i + f^/i , 

fti+1,2 = fti,2 + P,ifti,l + Ifi-l , 

fti+1,3 = fti,3 + Piifti,2 + Ifi-lf 11-1,1 , 

fli+l,fi ~^ f^tifli,fi—l ~^ T/i— l,At— 2 ; 

f^i+i,fi+i = PiJ.ffi,fj. + l^l-lf^l-l,^l-l ■ (34) 



The orthogonal polynomial and recurrence coefficients are recursively determined by ( |20D 
and (^3])-(0). The expansion coefficients for the optimized polynomial Pn{x) can be 
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obtained from 

= E /m- dxw{xff{x)x>'-' . (35) 



u=0 



The ingredients needed for this recursion are the basic integrals s^, defined in (|TT]) and 



A 



U = J dxw{xYf{x)x'' . (36) 

The recurrence scheme based on the coefficients of the orthogonal polynomials f^j^^ in 
(34) is not optimal for large orders ra, neither for arithmetics nor for storage requirements. 
A better scheme can be built up on the basis of the integrals 

r^,u = J dxw{xf^^{x)x'' , 

/i = 0, 1, . . . , n; u = fj,, n + 1, . . . ,2n — fj, . (37) 
The recurrence coefficients P^, ■y^-i can be expressed from 

= , P/. = r^,^+i + f^,lr^,|, (38) 

and eq. ( pOD as 

^ = -/.I - ^ , 7.-1 = ^ ■ (39) 

^fj-fj, '"/i— i,/-t— 1 

It follows from the definition that 

rX 

tqv = / dxw{xYx^ = Su , 

f-A 



riu = dxwix)\x''+' + fnx") = s,+i + fns, . (40) 
The recurrence relation (|19D for the orthogonal polynomials implies 

+ P^^rf,^ + 7^_ir^_i,;, . (41) 

This has to be supplemented by 

fu = -- (42) 
■So 



and by the first equation from (13^) 



U+i,i = U,i + P, . (43) 

Eqs. (^)-(|43|) define a complete recurrence scheme for determining the orthogonal poly- 
nomials $^(x). The moments s,y of the integration measure defined in (pUD serve as the 
basic input in this scheme. 
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The integrals bu in (pTsD, which are necessary for the expansion coefficients d^, in ([l2D , 
can also be calculated in a similar scheme built up on the integrals 



bf,u = / dxw{xff{x)^^{x)x'' 



fj, = 0,1, ... ,11] u = 0,1, ... ,n — fj, . (44) 



The relations corresponding to (^)-(|4l|) are now 

rX 



bou = J dxw{x) f{x)x'' = ty , 

rX 

bi^ = J dxw{xff{x){x''^'^ + fnx") = t^+i + futy , 

(45) 



The only difference compared to (|40|)-(^lD is that the moments of w{x) are now replaced 
by the ones of w{x)'^f{x). 

It is interesting to collect the quantities which have to be stored in order that the 
recurrence can be resumed. This is useful if after stopping the iterations, for some reason, 
the recurrence will be restarted. Let us assume that the quantities qi, (z/ = 0,...,n), 
bu iy = 0, . . . ,n), iy = 1, . . . ,n — 1) and ■ji, {y = 0, . . . ,n — 2) are already known 
and one wants to resume the recurrence in order to calculate these quantities for higher 
indices. For this it is enough to know the values of 

/n— 1,1 ; ^n— l,n— 1 ; 

^oln = ('"0,2n+l; 1^l,2nj ■ ■ ■ j 1^n,n+l) 5 
^oln = (^0,2n5 ''"l,2?i-l; • • • ; fn,n) ; 
Bo,ln = {bo,n, ^l,n-l5 • • • 5 &n,o) i 

-^O.-.n-l = (^0,n-lj bi^n-2, ■ ■ ■ , ^n-l,o) • (46) 

This shows that for maintaining a resumable recurrence it is enough to store a set of 
quantities linearly increasing in n. 

An interesting question is the increase of computational load as a function of the 
highest required order n. At the first sight this seems to go just like n^, which is surprising 
because, as eq. (^ shows, finding the minimum requires the inversion of an n^n matrix. 
However, numerical experience shows that the number of required digits for obtaining a 
precise result does also increase linearly with n. This is due to the linearly increasing 
logarithmic range of eigenvalues, as illustrated by (||). Using, for instance. Maple V for 
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the arbitrary precision arithmetic, the computation slows down by another factor going 
roughly as (but somewhat slower than) r? . Therefore, the total slowing down in n is 
proportional to n'^. For the same reason the storage requirements increase by r? . 



4 A convenient choice for TSMB 

In the TSMB algorithm for Monte Carlo simulations of fermionic theories, besides the 
simple function also the function x"°/-P(x) has to be approximated. Here P[x) is 
typically a lower order approximation to a;~". In this case, if one chooses to optimize the 
relative deviation, the basic integrals defined in (pf) and (^) are, respectively. 



A _ 2 2 + 



^ dx P{x)x''+'' . (47) 



It is obvious that, if the recurrence coefficients for the expansion of the polynomial P{x) 
in orthogonal polynomials are known, the recursion scheme can also be used for the 
evaluation of and t^,. 

Another observation is that the integrals in (|47|) can be simplifyed if, instead of choos- 
ing the weight factor w{xy = P(a;)^x^", one takes 



w[x 



|2 = P{x)x" , (4J 



which leads to 



A 



= ^ dxP{x)x 



A 



U = J dxx" . (49) 

Since P{x) is an approximation to the function f{x) = x~'^/P{x) is close to one and 
the difference between /(x)~^ and f{x)~^ is small. Therefore the least-squares optimized 
approximations with the weights w{x)^ = /(x)"^ and w(x)^ = f{x)~^ are also similar. It 
turns out that the second choice is, in fact, a little bit better because the largest deviation 
from typically occurs at the lower end of the interval x = e where P{x) is smaller 
than As a consequence, P{x)x" < 1 and P{x)x°' > P(x)^x^". This means that 

choosing the weight factor in (^) is emphasising more the lower end of the interval where 
P{x) as an approximation of is worst. 

In summary: least-squares optimization is a flexible and powerful tool which can serve 
as a basis for applying the two-step multi-bosonic algorithm for Monte Carlo simulations 
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of QCD and other similar theories. With the help of the recurrence scheme described in 
the previous section one can determine the necessary polynomial approximations to high 
enough orders. 
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Figure 1: The comparison of the relative deviation R{x) = xP{x) — 1 for the Chebyshev 
polynomial {Rc{x)) and the quadratically optimized polynomial {Ro{x)). In the lower 
part the two ends of the interval are zoomed. 
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a=l 




Figure 2: The (squared) deviation norm 5^ of the polynomial approximations oi x ^ as 
function of the order for different values of A/e. The asterisks show the e/A — > hmit. 
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Figure 3: The (squared) deviation norm 5^ of the polynomial approximations oi x "as 
function of the order at X/e— 10^ for different values of a. 
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